# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python

from __future__ import division
import scipy.weave as weave
from numpy import random, exp, sqrt, power, log, newaxis, zeros, \
        zeros_like, float32, int32, repeat, minimum, ones, mean, average, sum, \
        min, prod, nan, load

def creatematrices(p, w):
    """
    Creates rho, xi and delta 
    (as defined in the article)
    to compute the Rao-Blackwellised estimators.
    """
    rhomatrix = zeros((p+1, p+1))
    ximatrix = zeros((p+1, p+1))
    deltas = ones(p+1)
    code = \
    """
    float r;
    float temp;
    for(int i = 0; i < p + 1; i++){
        ximatrix(i, i) = 1.;
        for(int j = i + 1; j < p + 1; j++){
            if ((w(i) == 0) || w(j) > w(i) ){
                r = 1;
            }
            else{
                r = w(j) / w(i);
            }
            rhomatrix(i, j) = r;
            ximatrix(i, j) = ximatrix(i, (j-1)) * (1 - rhomatrix(i, j));
        }
        if (i > 0){
            temp = 0.;
            for (int j = 0; j < i ; j ++){
                temp += deltas(j) * rhomatrix(j, i) * ximatrix(j, i-1);
            }
            deltas(i) = temp;
        }
    }
    """
    weave.inline(code,["p", "w", "rhomatrix", "ximatrix", "deltas"], type_converters=weave.converters.blitz)
    return rhomatrix, ximatrix, deltas

def raoblackwell(k, p, omegaproposals, omegacurrent, weights):
    """
    Compute the RB2 weights leading
    to estimator hat tau 4 in the article.
    These weights correspond to phi in the article.
    """
    w = zeros(p+1)
    w[0] = omegacurrent
    w[1:(p+1)] = omegaproposals[:, k]
    rhomatrix, ximatrix, deltas = creatematrices(p, w)
    for i in xrange(0, p+1):
        weights[i, k] = deltas[i] * sum(ximatrix[i, i:(p+1)])
    weights[0, k] -= 1

